Local persistence in directed percolation 
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We reconsider the problem of local persistence in directed site percolation. We present improved 
estimates of the persistence exponent in all dimensions from 1 + 1 to 7+1, obtained by new 
algorithms and by improved implementations of existing ones. We verify the strong corrections 
to scaling for 2 + 1 and 3+1 dimensions found in previous analyses, but we show that scaling 
is much better satisfied for very large and very small dimensions. For d > 4 (d is the spatial 
dimension), the persistence exponent depends non-trivially on d, in qualitative agreement with the 
non-universal values calculated recently by Fuchs et al. (J. Stat. Mech.: Theor. Exp. P04015 
(2008)). These results are mainly based on efficient simulations of clusters evolving under the time 
reversed dynamics with a permanently active site and a particular survival condition discussed in 
Fuchs et al. These simulations suggest also a new critical exponent £ which describes the growth 
of these clusters conditioned on survival, and which turns out to be the same as the exponent, n + S 
in standard notation, of surviving clusters under the standard DP evolution. 



I. INTRODUCTION 

In general, persistence [H, H, E| in time-dependent crit- 
ical phenomena is the probability that some observable 
does not cross its long-time expectation until some finite 
time t. In particular, we will deal with the order pa- 
rameter, which in directed percolation is the density and 
satisfies (p) — > at the critical point. Local persistence 
P(r,t) is the probability that the local order parameter 
at position r does not cross this value up to time t. In 
directed percolation it is thus equal to zero when the site 
r was active in the initial configuration, while it is posi- 
tive and equal to the chance that it is not yet activated 
at t when it was inactive originally. In the following we 
shall only consider homogeneous systems in which case 
P(r,t) does not depend on r and will be written P(t). 

In general, local persistence decays according to a 
power law 



P(t) ~ r 6 



(1) 



where 9 is a new universal critical exponent, independent 
of the standard exponents. In particular, it is in general 
not related to the dynamical critical exponent z [l|, [2|, |3|, 
1- 

Previous studies of persistence (in the following we will 
only deal with local persistence) in the directed percola- 
tion (DP) universality class [j, [H, [|| have given some- 
what contradictory results. In particular, superuniver- 
sality was suggested in Q, i.e. the possibility that 9 is 
independent of dimension. This was later refuted in 
For d = 1, a connection with spreading of DP clusters 
in the presence of a wall wall was pointed out in [5[, al- 
though no direct relation between 9 and critical boundary 
exponents for DP could be established. 

In the following we shall again study this problem by 
means of numerical simulations. In spite of doubts con- 



cerning the universality of 9 Q, we restrict ourselves to 
site percolation on simple hyperbolic lattices. We use a 
new sampling strategy which is more time efficient than 
the standard strategy in low dimensions, and more space 
efficient in high dimension. In addition, we also made 
(for intermediate dimensions) simulations with a highly 
efficient implementation of the standard strategy. Fi- 
nally, since it was crucial to simulate exactly at the crit- 
ical point, we made new estimates of p c (the percolation 
threshold) in dimensions 1+2 to 1+7. One of our findings 
is that 9 not only is not superuniversal, but it depends 
also non-trivially on dimension above the critical dimen- 
sion for DP (which is d c = 4; in the following, we will 
always denote by d the spatial dimension). The latter 
had been predicted in Q, where it was also suggested 
that 9 might be not universal at all for d > d c . 

Our new sampling strategy suggests immediately a new 
critical exponent, called £ in the following. This expo- 
nent describes the growth of clusters under a modified 
evolution, where time is reversed, a permantent source of 
activity is added at the site at which persistence is mea- 
sured, and a rule is added that clusters die when they 
hit the baseline [1, [B[ . The growth of the mass of such 
clusters, conditioned on their survival, can be described 
by an exponent As pointed out by H. Hinrichsen after 
this paper was written Q , a heuristic argument relates £ 
to a combination of standard DP exponents which can be 
re- written, using hyperscaling and the standard notation 
for DP exponents [8[, as 



C = ry + S. 



(2) 



The r.h.s. is the exponent which describes the mass 
growth of surviving critical DP clusters under normal 
(forward time, single site seed) dynamics. 
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II. ALGORITHMS AND THE EXPONENT C 

In all simulations we used helical boundary conditions, 
i.e. a site is indexed by a single integer i e [0, . . ,N — 1], 
where TV is the size of the lattice. In all cases, N is 
chosen as a power of 2. Neighbors of i are indexed by 
i ± 1, i ± Li, . . . i ± id-i? all modulo AT. The integers 
are of order N 1 ^, but not necessarily equal to it (i.e. the 
lattice does not have to be an exact hypercube). 

The standard way to simulate persistence uses four 
data structures: 

• A bit (or, for implementation simplicity, byte) ar- 
ray of size TV, representing the activity pattern 
(sj = 1 if site i is active, Sj = otherwise) , which is 
replaced at each time step by a new empty array. It 
is used for checking whether a site to be activated 
is already active or not. 

• A list of indices of all active sites during the previ- 
ous time step. This list is used to activate neighbors 
in the next time step. 

• A similar list containing the active sites during the 
present time step, which originally is also empty 
and which will replace, after the time step is fin- 
ished, the previous list. 

• A bit array of size N representing the history of 
all sites: tj = 1, if site i had been active during 
any previous time step, and i, = otherwise. This 
array is updated continuously, but erased only at 
the start of a new run. 

The simulation proceeds then in the usual way, start- 
ing with the initial list of active sites (which can be a 
single site or a finite fraction of all sites) and simulating 
time step after time step, until either a maximal time is 
reached, all activity has died out, or until all sites have 
been activated at least once. For runs which start with 
a single active site, only a tiny fraction of the lattice will 
actually be activated when d is large. In that case a direct 
implementation of the bit arrays is wasteful of memory 
and is replaced by hashing. Details are discussed e.g. in 

While single-site starts are most efficient to locate the 
critical point and were used for this also in the 

present work, they cannot so easily be used for measuring 
persistence. For the latter, starts with a finite density of 
active sites (say p(0) — 1/2) are most straightforward, 
but they require enormous memory in high dimensions, 
since hashing is of no use. If we want to have no finite-size 
effects, we have to use lattices whose linear size is larger 
than t l / z , where t is the simulation time and z is the dy- 
namic critical exponent which for DP is z — V\\/v_\_ < 2. 
Even if we use multispin coding (i.e. 1 bit per site), 
we can only simulate rather short times on present-day 
workstations. 

A way out of this dilemma is to use the time-reversed 
process discussed in [1, H|. The persistence probability 



P(t) is defined as the probability that none of the space- 
time points (i, t') for fixed i and < t' < t is activated by 
any of the active sites on the initial hypersurface (i 1 , 0). 
This means that if site i is still persistent at time t, there 
cannot be any path of active sites (an analog argument 
holds also for bond percolation) connecting any of the 
sites (i,t') to any (i',0). But such paths, if they would 
exist, could also be followed in the opposite direction. 
Thus P{t) is also equal to the probability that the clus- 
ter of sites activated by "sources" on the line interval 
{(i,t'); —t<t'< 0} does not survive to time 0. 

This is implemented in the following recursive way (see 
Fig. 1). Let us denote by C t the cluster activated by 
the sites ("seed") {(i,t');-t < t' < 0}, and by H the 
hyperplane {(i',0);0 <= i' < N}. We start with t = 1, 
in which case Ct is just the site (i, — 1) plus all sites with 
t' > connected to it. If C t n Ho is not empty, i.e. if the 
seed activates at least one site in TCo, we discard Ct and 
start a new run. Otherwise, we add the site (i, — t — 1) 
to the seed and activate all sites connected to it which 
are not already in Ct- The new cluster, which will be 
at least as large as Ct but in general not much larger, is 
Ct+\ (or rather, due to the particularities of the lattices 
used in this paper, Ct+2] see the discussion at the end 
of this section). This is iterated until either the cluster 
intersects 7io or until t = t max is reached. 

It might seem at first that this is much more storage 
demanding than direct simulation, since we have to store 
now the entire space-time bit pattern, not only the spatial 
pattern at a fixed time. What saves us, however, is that 
the clusters are fractal with spatial fractal dimension < 2, 
whence hashing will be very efficient. 

While P(t) is the chance that Ct survives up to times 
> t, and thus 9 describes its survival probability, the 
new critical exponent C describes the growth of its mass. 
There are several ways to describe this growth. The first, 
and maybe most natural would be the ansatz 

M(t) = (#£> ~ t<S (3) 

where M(t) is the average mass of the still surviving clus- 
ters. Its main disadvantage is that contributions to M{t) 
from early times might not scale, but decrease only slowly 
in importance as t grows. 

Alternatively we can define critical exponents via the 
increase of integrated quantities like M(t), as this de- 
pends less on these non-scaling contributions. Thus in 
the following we will define £ via the average number 
m(t) of activated sites (except for the site at the tip of 
the cluster) during the step Ct-i — > Ct), 

m(t)~t c . (4) 

If there were no correlation between the cluster mass and 
its survival probability, we would have obviously 

Ci = C+ 1- (5) 

The same should still hold if there is exact scaling. Our 
numerics suggests that Eq. §5§ is correct, but with large 
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FIG. 1: (color online) Typical cluster of the time-reversed 
process, where a segment of a line i = 0, — t < t' < is the 
active seed. The growth is upward but in epochs, where each 
epoch corresponds to the downward extension of the seed by 
one site. The sites activated during the last epoch are given a 
different color. Notice that only every second site (arranged 
in a checkerboard pattern) is used for the process. Notice also 
that a straightforward growth of the cluster, without break- 
ing it up into epochs, would be much less efficient. Typical 
clusters would be very fat, before applying the conditioning of 
not cutting constant-i hyperplanes, and most clusters would 
not meet the condition. The main virtue of the present algo- 
rithm is that such clusters are eliminated already at a very 
early stage. 



fmite-i corrections. Appearant deviations from Eq.([5]) 
will dominate our error estimates. 

In addition to its mass we can also measure other char- 
acteristics of the cluster, such as its spatial extension. As 
usual in time dependent critical phenomena we can define 
z' as 

R{t)^t x / Z '. (6) 

The same definition can also be used for ordinary DP 
clusters, where z = v\\/v±_ [1] [2(J. We measured z' in 
d = 1 and d = 2 and found z' — z within rather small 
error bars (z' = 1.5803(9) for d = 1, z' = 1.770(9) for 
d = 2; the corresponding values of z are 1.5807 [ll| and 
1.767 [l2j]). Thus we conjecture that z' = z exactly. The 
fact that critical exponents describing geometric aspects 
are more robust than others is also well known from stan- 
dard critical phenomena, where entropic boundary expo- 
nents are different from bulk exponents, while correlation 
length exponents do not change near boundaries [l3j]. 

The relative efficiencies of direct simulations of half 
active lattices on the one hand, and the simulation of 



single time reversed clusters on the other hand, can be 
obtained by estimating the number of sites that have to 
be tested/activated to establish that one site persists up 
to time t. In the direct approach this is 

t 

nStod«Ep(*)/i'(*)~* fl+1 " a . (7) 

{'=0 

where 8 is the exponent for the decay of the density of 
active sites. For the time reversed cluster growth it is 

«eT~* C+1 - (8) 

The ratio is 

reverse / direct +(+5-8 (q\ 

"tested / "tested 1 l y J 

If the exponent in this formula is negative, the time re- 
versed cluster growth has less time complexity than the 
direct simulation. 

Before leaving this section, let us make three remarks: 

• For the estimates of p c we used the variance re- 
duction method described in [3, [H[. This gave 
particularly big improvements for d > 4. 

• As noticed already in 0, @, the regions left and 
right of the site at which persistence is measured 
are decoupled and can be treated independently. 
Thus the sites i' > and i' < in Fig. [T] can 
be simulated independently. Both sides contribute 
the same to M(t) and to logP(i), so decay is much 
slowed down when only one half of the cluster is 
simulated and accuracy is substantially improved. 

• For any d we use lattices where each bond changes 
only one of the spatial coordinates and, of course, 
time (more precisely, an active site can activate 
nearest neighbors in a simple hypercubic lattice at 
exactly one unit later time). Thus the lattices sep- 
arate naturally into two checkerboard type sublat- 
tices. If activation is at the start restricted to one of 
them, it stays on it forever. Thus we can, without 
loss of generality, start with active sites restricted 
to one of the sublattices (for d = 1 this is clearly 
seen in Fig. [TJ. Again this divides the exponent 
9 by a factor 2 as compared to simulations where 
both sublattices are active, and makes simulations 
more easy. The values of 9 quoted below refer to 
activation restricted to one sublatticc. 



III. RESULTS 

Our main results are summarized in Table 1. Except 
for d = 1, where extremely precise estimates of p c are 
available from series expansions [llT | , we first made stan- 
dard spreading simulations [Toj where we measured the 
mass, survival probability, and r.m.s radius of clusters 
grown from single point seeds. For d = 4 we took into 
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deviations from scaling noticed already in [J]. A priori, 
the strongest deviations would have been expected for 
d = 4 since this is the upper critical dimension for DP, 
and we must expect logarithmic corrections there as for 
all other observables. It is completely unclear why the 
corrections to scaling are largest for d = 3 instead, and 
why they also seem to decrease slowest with t (they seem 
better compatible with logarithmic corrections than with 
power-law terms). Indeed, the data for d = 3 taken by 
themselves might suggest a different relation between 9 
and t than a power law, e.g. a stretched exponential. 
Although extremely unlikely from a theoretical point of 
view, we cannot really exclude this possibility. For d > 4 
the corrections to scaling are still large, but exponents 
defined for them via 

P(t) ~ t- s (l + (10) 

seem to increase. Thus, in spite of the visible curvature 
of all curves in Fig. 2, determination of 9 becomes more 
reliable for larger d, as soon as d > 5. 

Values of 9, with subjective error bars dominated by 
the uncertainties of parameterizing the scaling correc- 
tions, are shown in Table 1. While our estimates of 9 are 
typically larger than those of @ , they agree with those of 
[J] with one minor exception: While 9 > 1.62 is quoted 
in [|[ for the contact process in d — 2 (and extrapolating 
the corresponding curve in Fig. 3 of their paper would 
suggest indeed a value substantially larger than 1.62), we 
obtain 1.611(1) for d = 2. We believe that this is most 
likely due to an inaccurate value of the critical rate A c 
used in 

Table 1 suggests that 9 — ► 1 in the limit d — > oo, as pre- 
dicted analytically in Q . But the agreement is not quan- 
titative, as our values are much larger than those given 
by Eq. (14) of 4]. It might be that this is because the 
latter was derived for bond percolation, but this seems 
unlikely in view of the good agreement between expo- 
nents obtained for (site) DP and for the contact process. 
This agreement strongly suggests universality. 

Log-log plots of the size M(t) of inverse dynamics clus- 
ters Ct against t are shown in Fig. [3) To stress that £i = 2 
within statistical errors for d > 4, and to reduce the plot- 
ting range on the y-axis, we plotted actually M(t)/t 2 . 
Notice also that points on the central line (which belong 
to Ct trivially) are excluded from M(t). We see strong 
deviations from pure power laws for all dimensions 2, 
but the data are in all cases compatible with asymptotic 
power laws. Similarly large deviations from pure power 
laws are also seen in log-log plots of m{t) versus t (sec 
Fig. [4j, but the correction to scaling exponents A, de- 
fined via Eq. (|10[) . are systematically larger there. Thus 
extracting asymptotic power laws from Fig. [¥] is easier 
than from Fig. [31 in spite of larger statistical fluctua- 
tions. The asymptotic exponents extracted from Figs. [3] 
and|3]should be related by Eq.©. For all dimensions this 
relation is roughly fulfilled, with the resulting estimates 
quoted in Table 1. The errors quoted there are mostly 



TABLE I: Main results. All values of p c are new estimates, 
except the value for d = 1 which is from [ll|] . £ max is the maxi- 
mum time over which clusters of the time-reversed process are 
grown, and N c i is the number of clusters grown in this way, 
including those which died before reaching £ max . For d = 2 
also a large number of direct simulations were made, typically 
up to the same £ max , and their results are also taken into ac- 
count in the estimate of 9. The question mark for 9(d = 3) 
indicates that the data are also compatible with some other 
relation (e.g. a stretched exponential) instead of a power law. 
The last column uses values from Hi. 
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FIG. 2: (color online) Estimates of t 15 P(t) against t, for d = 1 
(bottom) to d = 7 (top). Notice the very large corrections to 
scaling for d = 3 which are indeed larger than those for d = 4 
where we expect logarithmic corrections. 



account the logarithmic corrections calculated in [1 61 ] in 
the same way as in [3]. Otherwise, we estimated the 
value of p c by demanding that the total number of all 
active sites is a pure power law for large t, up to possible 
corrections to scaling. Virtual lattice sizes were in each 
case 2 64 , which allowed for t values at least twice as long 
as those used for measuring persistence. The results are 
given in the second column of Table 1. 

Results for P(t) are shown in Fig. 2. The curves for d = 
2 to 7 show the raw data (multiplied by t 1 - 5 ), while the 
curve for d = 1 shows the square of the actually measured 
survival probability (see previous section). The data for 
d = 1 show by far the smallest corrections to scaling - 
they are well represented by a power law for t > 100. In 
contrast, for d = 2 and d = 3 we see the extremely strong 
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FIG. 3: (color online) Estimates of M(t)/t 2 against t, for 
d — 1 (bottom) to d = 7 (top). Here M(t) is the number of 
active space-time sites (excluding the center line) in inverse 
dynamics clusters Ct- 
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FIG. 4: (color online) Estimates of m(t)/t against t, for d — 1 
(bottom) to d — 7 (top). Here m{t) is the number of active 
sites (excluding the center line which is always active) added 
to the inverse dynamics clusters Ct during the step [t— 2) — > t. 
To reduce statistical fluctuations, data have been binned with 
At/t » 0.03. 



reflecting the discrepancies between the estimates based 
on Figs. [3] and |H 

After submission of this paper, it was pointed out to me 
by Haye Hinrichsen [7| that C can be related to standard 
exponents. If we assume that the density in the time 
reversed cluster decays radially, from a value near unity, 
as r~ zS with a radial cutoff at r> ~ t 1 ^, then 



( = d/z-5 . 



(11) 



Using the hyperscaling relation d/z = 26 + /i for d < 4, 
where r\ describes the average number of sctive sites in 
single-site seeded DP clusters [8|, this can be written as 



C = V + 5 



(12) 



While Eq. (HI]) can only hold for d < 4 (for d > 4 the 
assumption of a power law radial density decay from a 
value of order unity must fail), Eq. (fT2")) seems to hold 
for all dimensions. To demonstrate this, we quote in the 
last column of Table 1 the values of ij + 5 from [8| . In all 
cases, a greement holds to better than two standard de- 
viations. Notice that Eq. (fP2"|) is very surprising. It says 
that two cluster growth mechanisms which seem to have 
no close relationship, and for which the cluster survival 
probabilities scale with different critical exponents, have 
nevertheless the same growth exponents after condition- 
ing on survival. 

Finally, let us discuss the relative time complexities of 
direct and time reversed (single cluster) simulations us- 
ing Eq.([9]) and the critical exponents given in Table 1 
and in [8| . For d = 1 we find that the time reversed sim- 
ulation is much more efficient, by a factor ~ t . Thus 
one gains several orders or magnitude in CPU time. The 
same is qualitatively true for d — 2, although the differ- 
ence is much less, only a factor ~ t° . For d — 3 both 
methods have the same time complexity, and for d > 3 
the ranking is reversed - quite substantially (t~ - 9 ) for 
d = 7. Thus direct simulations would be much faster in 
high dimensions, if we could solve the space complexity 
(memory limitation) problem. The other advantage of 
time reversed simulations is of course that they offer the 
unique possibility to measure the exponent 



IV. DISCUSSION 

We presented in this paper numerical estimates for the 
local persistency exponent 9 in DP. These estimates are 
more precise than previous ones for spatial dimensions 
d < 4, and are the first published ones for d > 4. The 
latter were made possible by a combination of hashing 
(virtual memory) and a novel way to simulate the growth 
of the time reversed clusters discussed in The latter 

algorithm suggested also very naturally another new crit- 
ical exponent, £. It seems that 9 and Q are not trivially 
related, but that £, which describes the cluster growth 
of an artifical time-reversed evolution with modified dy- 
namics, is the same as the usual DP exponent describing 
mass growth of active clusters. While this is supported 
both by an indirect heuristic argument and by numerics, 
it seems very surprising and not well understood intu- 
itively. 

We verified the very large corrections to scaling seen in 
previous analyses local persistency in DP, in particular 
for 2 and 3 spatial dimensions. Indeed, these corrections 
seem to be larger in d = 3 than in d — 4, where we 
would have expected a priori that logarithmic corrections 
should give the strongest deviations from a pure power 
law. Although we see no theoretically appealing scenario 
for it, the data alone would suggest that the asymptotic 
behavior of P(t) in three dimensions is not described by 
a power law at all, but rather by a stretched exponential 
or something similar. 
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In the present paper we dealt only with directed per- 
colation, but very similar problems can be studied also 
for ordinary (undirected) percolation. The situation is 
maybe best described in two dimensions (analogous to 1 
spatial plus 1 temporal dimension in DP). Let us con- 
sider a wedge cut out from a square lattice, with edges 
W\ and W% meeting in point P and extending to infinity 
away from P. The angle between W\ and W2 is a. Let us 
furthermore consider intervals 1\ and I2 on W\ resp. W2, 
with lengths l\ and l^. We can then ask for the prob- 
ability P(Ii, I2, a) that there is no path from any point 
on 1\ to any point on I2. This is complementary to the 
probability that there is a path from I\ to I2 ■ Probabili- 



ties of the latter type (called 'crossing probabilities') were 
studied extensively [13, HH, Eil , but it seems that neither 
P(Ii, I2, a) nor its generalizations to higher dimensions 
were studied previously. 
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